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0.1  Objectives 

The  United  Technologies  Research  Center  (UTRC),  a  Division  of  UTC,  was  funded  by 
AFOSR  Computational  Mathematics  Program  to  investigate  a  class  of  singularities  that  arise 
in  mathematical  models  describing  interface  dynamics  in  multiphase  flow.  The  main  focus 
was  the  process  of  liquid  ligament  pinch-off,  the  understanding  of  which  is  crucial  to  the 
first-principle  modeling  of  fuel  spray  atomization.  The  singularities  arise  from  the  continuum 
Navier-Stokes  formulation  of  two-phase  flow  as  liquid-gas  interfaces  approach  each  other 
towards  pinch-off.  Physically,  the  continuum  concept  of  surface  tension  at  an  infinitesimally 
thin  interface  does  not  hold  when  molecular  fluctuations  begin  to  be  comparable  to  inertial 
and  viscous  forces.  The  question  to  be  addressed  is  whether  there  is  a  unique  continuation 
(through  singularity)  of  the  continuum  solution  obtained  from  the  model  discretization 
implemented  in  multiphase  algorithms.  From  the  computational  standpoint,  the  crucial  point 
is  whether  initial  conditions  after  the  singularity  depend  on  the  details  of  the  numerical 
discretizations  leading  to  pinch-off. 

The  objective  of  the  project  is  to  develop  a  physically/mathematically  consistent  multi-scale 
computational  approach  to  continue  the  solution  through  the  formal  singularity  described 
above.  The  multi-scale  approach  consists  of  a  macroscopic,  time-resolved  interface-capturing 
simulation  of  two-phase  flow;  a  singularity-free  mesoscopic  simulation  that  bridges  atomic 
and  continuum  scales;  and  a  physics-based  closure  model  derived  from  the  mesoscopic 
analyses  to  be  embedded  into  the  macroscopic  continuum  to  bypass  the  singularity  issue. 


0.2  Summary  of  accomplishments 

The  overall  objective  of  developing  a  mathematically/physically  consistent  multiscale  and 
multiphase  numerical  approach  to  continue  simulations  through  singularities  due  to 
topological  changes  of  the  interface  was  accomplished.  The  simulation  framework  developed 
under  this  project  will  be  applied  to  future  high-fidelity  simulations  of  engine  spray 
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atomization  to  accurately  predict  fuel  droplet  statistics  and  guide  engine  design,  and  it  is 
expected  to  have  a  substantial  impact  in  the  jet-engine  industry. 

Towards  the  first  objective  of  documenting  and  explaining  the  inconsistencies  of  continuum 
direct  numerical  simulations  of  non-disperse  two-phase  flow,  it  was  successfully 
demonstrated  that  the  mathematical  singularity  of  assuming  an  infinitesimally  thin  interface 
near  the  pinch-off  of  a  liquid  jet  leads  to  divergence  of  simulation  in  terms  of  interface 
position  and  flow  pressure.  An  advanced  interface-capturing  multiphase  solver  developed 
previously  by  UTRC  was  used  to  simulate  a  canonical  case  of  drop  pinch-off  under  capillary 
and  straining  effects.  After  the  lack  of  convergence  with  grid  refinement  was  demonstrated, 
several  algorithms  were  tested  to  bypass  the  pinch-off  singularity  by  numerically  forcing 
breakup  at  an  assigned  threshold.  The  results  are  documented  in  Subsection  1.1. 

To  achieve  a  better  understanding  of  interface  fluctuations,  the  project  next  was  focused  on 
developing  a  mesoscopic  simulation  approach  that  bridges  molecular  and  continuum  scales. 
The  development  was  based  on  the  general  framework  of  dissipative  particle  dynamics 
(DPD).  Several  accomplishments  are  worth  mentioning.  A  many-body  dissipative  particle 
dynamics  (MDPD)  code,  obtained  by  combining  short-range  repulsive  and  long-range 
attractive  forces,  was  applied  to  vapor/liquid  simulations.  Based  on  the  radial  distribution  of  the 
virial  pressure  in  a  drop  at  equilibrium,  a  systematic  study  was  carried  out  to  characterize  the 
sensitivity  of  the  surface  tension  coefficient  with  respect  to  the  inter-particle  interaction 
parameters:  particularly,  the  approximately  cubic  dependence  of  the  surface  tension  coefficient 
on  the  bulk  density  of  the  fluid  was  evidenced  for  the  first  time.  In  capillary  flow,  MDPD 
solutions  were  shown  to  satisfy  the  Plateau-Rayleigh  condition  on  the  wavelength  of  an  axial 
disturbance  leading  to  the  pinch-off  of  a  cylindrical  liquid  thread;  correctly,  no  pinch-off 
occurred  below  the  cutoff  wavelength.  The  results  are  documented  in  subsection  1.2. 

The  objective  of  developing  a  coupling  between  the  particle  (P)  and  continuum  (C) 
formulation  was  limited  to  the  one-way  coupling  of  the  continuum  simulation  providing 
time-dependent  boundary  conditions  to  the  particle  simulation.  This  required  the 
development  of  a  Non-Periodic  Boundary  Condition  (NPBC)  algorithm  for  two-phase  flow 
that  can  preserve  a  liquid  gas  interface  through  a  non-periodic  boundary.  Several  non¬ 
periodic  algorithms  are  described  in  the  literature  for  bulk  flow  fluid  dynamics,  but  none,  to 
the  Pis’  knowledge,  for  free  surface  flow.  The  full  two-way  coupling  was  found  to  be 
irrelevant  to  the  overall  objective  of  developing  a  closure  model  to  allow  continuation 
through  singularity.  Particularly,  it  was  found  that  the  objective  of  two-phase  flow 
applications  only  require  small-scale  information  to  be  fed  into  the  large-scale  simulation, 
and  not  vice-versa.  Thus  only  one-way  coupling  is  required.  The  funding  was  then  focused 
on  the  development  of  a  physics-based  closure  model  and  a  generic  implementation  of  the 
model  in  the  continuum  direct  numerical  simulation. 

To  develop  a  continuum  closure  model  based  on  detailed  physics,  scaling  analysis  was 
performed  on  the  mesoscopic  simulation  results  from  MDPD.  The  analysis  illustrates  the 
cascade  of  fluid  dynamics  behaviors  from  potential  to  inertial-viscous  to  stochastic  flow,  and 
the  dynamics  of  a  liquid  jet  minimum  radius  is  consistent  with  the  power  law  scaling 
predictions  from  asymptotic  analysis.  One  major  accomplishment  is  the  key  finding  that  the 
scaling  transition  is  dependent  on  fluid  properties  only  and  independent  of  external  flow 
conditions.  The  transition  point  can  therefore  be  used  as  the  threshold  to  enforce  numerical 
breakup,  leading  to  physically  accurate  and  grid  converged  solutions.  Because  of  the 
condition-independent  nature  of  the  threshold,  the  same  threshold  for  breakup  can  be  applied 
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for  any  pinch-off  processes  occurring  in  the  simulation  domain,  independently  from  their 
local  strain  environment.  Thus,  the  scaling  transition  point  becomes  the  only  parameter  for 
the  closure  model.  This  model  was  implemented  in  the  continuum  two-phase  flow  solver  in  a 
completely  grid-independent  fashion.  The  implementation  was  demonstrated  to  converge 
under  the  physical  breakup  threshold  and  thus  provide  physically,  mathematically  and 
computationally  consistent  results.  These  results  are  documented  in  subsection  1.3  and  1.4. 


0.3  Organization  of  the  report 


A  comprehensive  set  of  results  is  available  in  published  papers  and  from  the  proceedings  of 
conferences  attended  by  the  Pi’s.  In  this  report,  we  provide  more  details  for  the  key  results 
obtained  during  the  project.  References  to  journal  papers,  conference  papers  and 
presentations  are  provided.  A  list  of  references  is  available  in  the  last  chapter  of  the  report. 


3 

Distribution  Statement  A:  Approved  for  public  release;  distribution  is  unlimited 


Contents 


0.1  Objectives  .  1 

0.2  Summary  of  accomplishments .  1 

0.3  Organization  of  the  report .  3 

1  Summary  of  research  results  5 

1 . 1  Demonstration  singularity  in  continuum  two-phase  simulations .  6 

1.1.1  Continuum  two-phase  approach .  6 

1.1.2  Singularity  demonstration  in  continuum  simulations .  6 

1.1.3  Bypassing  singularity  by  threshold  breakup:  numerical  tests .  9 

1 .2  Development  of  a  singularity-free  meso-scale  simulation  approach .  12 

1.2.1  B  ackground  ofMDPD .  12 

1.2.2  Numerical  details  of  MDPD .  13 

1.2.3  Application  to  liquid-vapor  interface  dynamics .  15 

1 .2.4  Non-periodic  boundary  conditions .  19 

1.3  Scaling  analysis  of  pinch-off  using  meso-scale  simulations .  22 

1.3.1  Theoretical  background  of  pinch-off  scaling .  22 

1 .3.2  MDPD  scaling  analysis  for  capillary  pinch-off .  23 

1.3.3  MDPD  scaling  analysis  for  pinch-off  in  straining  flow .  24 

1 .4  Development  and  demonstration  of  continuum-scale  closure  model .  27 

1.4.1  Closure  model  development  from  scaling  analysis .  27 

1.4.2  Generic  implementation  of  the  closure  model .  27 

1.4.3  Demonstration  of  simulation  convergence  using  the  closure  model ....  29 

2  Personnel  supported  31 

3  Publications  and  presentations  32 


4 

Distribution  Statement  A:  Approved  for  public  release;  distribution  is  unlimited 


Chapter  1 

Summary  of  research  results 


While  multiphase  flow  involving  complex  topological  changes  occurs  in  a  wide  range  of  Air 
Force  applications,  observation  of  the  multi-scale  interfacial  dynamics  in  a  laboratory  is 
limited  by  the  difficulty  in  controlling  the  boundary  conditions  of  the  experiment  and  by 
constraints  on  instrument  time  response  versus  available  field  of  view.  The  development  of  a 
mathematically  correct  and  computationally  feasible  approach  is  a  key  step  to  achieve  a  deep 
understanding  of  the  process. 

This  project  aims  to  address  a  fundamental  limitation  in  the  first-principle  simulations  of 
multiphase  flow.  Asymptotic  analysis  indicates  that  the  Navier-Stokes  (NS)  equations  for 
multiphase  flow  lead  to  interface  singularities.  In  this  report,  the  limitation  of  simulations 
based  on  the  discretization  of  multiphase  NS  equations  is  first  demonstrated  using  a 
macroscopic  Coupled  Level  Set  and  Volume  Of  Fluid  (CLSVOF)  interface-capturing 
approach  (subsection  1.1).  To  identify  when  microscopic  dynamics  affects  macroscopic 
physics,  a  computer  code  based  on  Many-body  Dissipative  Particle  Dynamics  (MDPD)  was 
developed  that  bridges  molecular  and  continuum  scales  (subsection  1.2).  This  effort  enabled 
the  creation  of  a  “virtual  test  facility”  for  small-scale  multiphase  flow,  whose  results  do  not 
depend  on  prior  calibration.  The  “facility”  allows  the  derivation  of  detailed  asymptotic  results 
for  pinch-off  (subsection  1.3)  and  the  development  of  a  closure  model  embedded  in  the 
macroscopic  simulation  to  achieve  physically/mathematically  consistent  results  in  the 
presence  of  liquid  pinch-off  (subsection  1.4).  The  main  elements  of  this  research  are 
summarized  in  the  sketch  below. 


Macroscopic  Mesoscopic 


Figure  1.  Summary  of  research  elements  for  the  investigation  of  computational 
continuation  through  multiphase  singularity. 
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1.1  Demonstration  of  singularity  in  continuum  two-phase 
simulations 


1.1.1  Continuum  two-phase  approach 

In  the  continuum  approach  adopted  for  this  project,  the  CLSVOF  method  [1]  is  used  to 
directly  capture  the  evolution  of  a  liquid-gas  interface.  Briefly,  the  Navier-Stokes  equations 
for  incompressible  flow  of  two  immiscible  fluids  are  augmented  with  a  smooth  level  set 
function  <j>,  whose  zero  level  represents  the  time-evolving  interface.  In  addition  to  the 
evolution  equation  for  (f),  the  transport  equation  for  the  cell  liquid  volume  fraction  (the 
volume-of-fluid  function,  F)  is  also  solved.  The  interface  normals  used  in  the  VOF 
reconstruction  step  are  determined  from  the  level  set  function.  The  volume  fractions  and  the 
normals  are  then  used  to  construct  a  volume  preserving  distance  function  (j).  Essentially, 
volume  is  preserved  by  implementing  a  “local”  mass  correction  at  every  iteration.  Second- 
order  accurate  curvature  is  calculated  from  F  by  the  method  of  height  fraction. 

The  solver  is  implemented  under  the  framework  of  adaptive  mesh  refinement  (AMR)  to 
reduce  the  computational  cost  of  resolving  multi-scale  interface  dynamics.  A  detailed 
description  of  the  block-structured  AMR  can  be  found  elsewhere  [2],  Cells  that  are  crossed 
by  the  liquid-gas  interface  are  tagged  for  refinement.  Starting  from  the  base  level,  boxes  (with 
a  minimum  size  of,  say,  323  cells)  are  combined  to  cover  all  the  tagged  cells  within  assigned 
coverage  efficiency.  This  set  of  blocks  with  the  same  grid  spacing  forms  level  1.  This  level  is 
in  turn  tagged  for  refinement  at  the  interface,  and  the  process  is  repeated  until  the  required 
grid  resolution  is  achieved.  The  block-structured  AMR  framework  allows  easy 
implementation  of  special  finite  differencing  schemes  for  multiphase  flow  on  Cartesian  grids. 

1.1.2  Singularity  demonstration  in  continuum  simulations 

We  consider  the  pinch-off  of  a  dumbbell-shaped  drop  due  to  either  capillary  effects  or 
elongation  effects  in  a  straining  flow.  The  simulation  convergence  is  tested  by  a  grid- 
refinement  study  using  AMR.  The  lack  of  convergence  in  the  simulations  shown  next  points 
to  the  singularity  due  to  the  topology  change  of  the  drop.  The  axisymmetric  water  drop  is 
immersed  in  an  extensional  flow,  as  shown  in  Fig.  2,  with  velocity  field 


ur(r,x)  =  - yr  12 
ux(r,x)  =  yx 


where  f  is  the  straining  rate.  The  center  of  the  drop  overlaps  with  the  stagnation  point  of  the 
extensional  flow.  Its  initial  shape  is  prescribed  as 


r(x)  = 


0.5 D  -  £XOS(2 7TX  /  X) 
^(().5D)2  -(\x\-0.5Af 


|jc|  <  0.5/1 

0.5/1  <|x|  <0.5/1  +  0.50 


In  the  absence  of  the  external  straining  flow,  the  drop  may  experience  breakup  driven  by 
capillary  effects  or  recover  to  spherical  shape,  depending  on  the  ratio  A/D  being  larger  or 
smaller  than  n  [27],  With  extensional  straining,  the  drop  can  be  elongated  and  finally  broken 
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by  the  external  flow  if  the  strain  rate  is  large  enough.  Here  we  select  two  different  values  of 
wavelength,  under  which  the  drop  breaks  up  by  capillary  effects  QJD=3.75>n)  and  by 
elongation  at  high  strain  rate  ’k/D=2.5<n ,  y  =  1 03  s'1 ),  respectively. 


x 


A 


Figure  2.  The  dumbbell-shaped  water  drop  placed  in  the  center  of  an  extensional 
straining  air  flow. 

In  the  following  tests,  the  computational  domain  is  set  to  be  15  cm  x  60  cm  and  covered  with 
a  uniform  base  grid  of  64  x  256.  Refinement  was  performed  by  successively  adding 
refinement  levels  near  the  water- air  interface,  yielding  effective  grid  resolutions  of  128  x  512 
at  one  level  and  256  x  1024  at  two  levels. 

The  locations  of  the  water-air  interface  prior  to,  near,  and  post  pinch-off  are  shown  in  Fig.  3 
for  the  capillary  case  ( y  —  0 ),  in  (a),  and  for  the  straining  case  (y  =  103  s'1 ),  in  (b).  Results 
from  different  refinement  levels  are  shown  together  in  the  same  plots.  It  can  be  observed  that 
the  interface  locations  converge  prior  to  pinch-off  for  both  cases.  Near  pinch-off, 
discrepancies  between  different  refinement  levels  are  observed,  and  even  larger  discrepancies 
are  found  post  pinch-off.  Comparing  the  two  cases  in  (a)  and  (b),  we  find  that  the  breakup 
under  the  straining  flow  is  more  sensitive  to  grid  discretizations  for  this  set  of  conditions.  The 
interface  disconnects  due  to  numerical  interface  capturing  issues  when  the  neck  region 
becomes  under-resolved. 
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Figure  3.  The  interface  location  predicted  using  different  resolutions  for  the  capillary 
case  (a)  and  extensional  case  (b).  Corresponding  maximum  pressure  at  the  pinch-off 
location  for  the  capillary  case  in  (c)  and  the  extensional  case  in  (d). 

The  corresponding  maximum  pressure  at  the  pinch-off  location  is  plotted  in  Fig.  3(c)  for  the 
capillary  case  and  in  Fig.  3(d)  for  the  straining  case.  Both  sets  of  curves  show  that  the 
maximum  local  pressure,  as  a  function  of  time  reaches  a  peak  at  pinch-off.  At  different  grid 
resolutions,  the  peak  values  differ,  increasing  with  decreasing  grid  size  because  of  smaller 
ligament  width  that  can  be  resolved  by  the  grid  and  therefore  higher  interface  curvatures.  It 
can  be  extrapolated  that  the  pressure  will  not  converge  with  further  increase  of  resolution 
because  the  theoretical  value  of  pressure  approaches  infinity  at  pinch-off.  Because  of  the  lack 
of  convergence,  the  prediction  of  pinch-off  time  becomes  highly  sensitive  to  grid  resolution. 
As  shown  in  Fig.3  (c)  and  (d),  the  peak  pressure  shows  up  at  different  times  for  different  grid 
resolutions  and  the  discrepancies  are  further  amplified  in  the  extensional  case  (d)  compared 
to  the  capillary  case  (c). 

The  above  simulation  results  tell  that  a  grid-converged  simulation  of  pre-pinch-off  liquid-gas 
interface  diverges  when  the  interface  approaches  pinch-off  point.  This  result  is  consistent 
with  the  Laplace  law  predicting  that  the  capillary  pressure  inside  the  liquid  ligament  depends 
linearly  on  the  liquid  surface  curvature.  Numerically,  the  high-pressure  build  up  inside  the 
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ligament  leads  makes  it  increasingly  difficult  for  the  solver  to  provide  an  accurate  flow 
solution.  A  closure  model  is  therefore  necessary  to  continue  simulating  the  post  pinch-off 
portion  and  bypass  such  computationally  challenging  point. 


1.1.3  Bypassing  singularity  by  threshold  breakup:  numerical  tests 


Before  introducing  the  physics-based  closure  model  in  subsection  1 .4,  we  verify  our  findings 
on  singularity  by  numerical  forcing  breakup  at  a  prescribed  threshold  to  bypass  the 
singularity  and  achieve  simulation  convergence. 


P->o/Rc  grad(cp)->0  Multiple  band  t-^c 

cell  tagging 

(a)  (b)  (c)  (d) 

Figure  4.  Different  numerical  approaches  to  introduce  threshold  breakup. 

For  the  time  being,  the  threshold  is  set  for  demonstration  purposes  and  we  do  not  apply  any 
physical  argument  to  determine  its  value.  The  most  critical  step  in  forcing  breakup  is  to  find 
the  instant  when  the  ligament  reaches  the  prescribed  threshold.  Once  the  simulation  reach  this 
threshold  length  l  f  the  state  at  the  center  of  the  ligament  can  be  forced  to  change  from  liquid 

to  gas  by  modifying  the  level  set  and  volume  of  fluid  functions.  As  shown  in  the  schematics 
in  Fig.  4,  there  are  several  different  ways  of  defining  t  f  and  correspondingly  different 

numerical  approaches  to  force  breakup. 

The  first  approach,  shown  in  frame  (a),  uses  a  value  of  critical  pressure.  For  a  thin  cylindrical 
ligament,  the  pressure  difference  between  liquid  inside  and  gas  outside  is  ZlP  =  cr/r.  When 
the  ligament  size  reaches  the  threshold  length  l  f  ,  the  pressure  inside  reaches  a  critical  value 

Pc  =y/£f  +Pg  -  assuming  the  pressure  in  the  gas  phase  Pg  is  constant.  For  each  step  in  our 
simulation,  we  search  for  a  location  that  has  a  global  maximum  pressure  value.  If  the 
maximum  pressure  reaches  the  critical  pressure  P  ,  we  force  breakup  to  occur  at  this 
location. 

The  second  approach  in  frame  (b)  and  the  third  approach  in  frame  (c)  are  based  on  the  level 
set  geometric  description  of  the  interface.  In  the  level  set  description,  a  distance  function  (f>  is 
used  to  define  a  field  that  has  the  interface  embedded  as  the  zero  distance  contours. 
Therefore,  the  distance  to  the  interface  is  explicitly  available  at  any  location.  A  narrow-band 
with  size  i  f  /2  can  be  defined  so  that  the  distance  from  all  the  points  within  this  band  to  the 

interface  is  less  than  f//2.  The  thinning  process  of  a  ligament  can  be  viewed  as  the  l  f  / 2- 
size  band  contracts  and  moves  inward.  When  the  band  merge  with  itself  at  a  certain  location, 
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the  size  of  the  ligament  is  decreased  to  l  f  and  the  ligament  is  numerically  forced  to  breakup. 

We  found  two  numerical  approaches  to  identify  the  location  and  instant  for  band  merging. 
One  as  shown  in  Fig.  4(b)  is  to  find  a  cell  such  that  the  module  of  the  gradient  of  the  level  set 
function,  I  V^l,  significantly  deviates  from  1  (For  a  level  set  function  defined  by  distance  to  a 
surface  with  radius  of  curvature  larger  than  f//2,  |V^|  =1  always  holds).  Another  approach, 
shown  in  Fig.  4(c)  is  to  find  a  cell  that  is  tagged  as  a  band  cell  for  multiple  times. 

The  fourth  approach  shown  in  Fig.  4(d)  is  the  direct  prescription  of  the  pinch-off  instant 
irrespective  of  the  ligament  state.  This  approach,  requiring  a-prior  knowledge  of  the  pinch-off 
location  and  therefore  not  applicable  to  a  generic  configuration,  is  expected  to  provide 
reference  results  with  the  least  amount  of  numerical  discrepancy  at  different  resolutions. 


(c)  (d) 

Figure  5.  The  interface  location  predicted  at  different  resolutions  using  different 
threshold  breakup  approaches  as  in  Fig.  4. 

The  interface  location  and  maximum  pressure  for  the  straining  case  (y  =  \  03  s'1 )  using  the 
four  different  threshold  breakup  approaches  are  shown  in  Fig.  5  and  Fig.  6  respectively. 
Results  from  different  refinement  levels  are  shown  together  in  the  same  plots.  Compared  to 
Fig.  3  (b)  and  (d),  all  the  threshold  approaches  predict  much  less  discrepancy  between 
simulations  at  different  resolutions  near  and  post  pinch-off.  Due  to  the  numerical 
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approximation  involved  in  computing  the  threshold  pressure,  the  level  set  gradient  and  level 
set  band  width  near  the  pinch-off  location,  different  approaches  reflected  different  type  of 
sensitivities  to  the  grid  resolution.  Even  if  using  the  same  ligament  width  threshold,  the  level 
set  gradient  approach  in  frame  (b)  and  the  multi-band  cell  tagging  approach  in  (c)  predict 
different  convergence  behaviors.  As  expected,  the  pinch-off  time  prescription  in  (d)  predicts 
the  least  amount  of  discrepancy  between  different  resolutions. 


(a)  (b) 


(c)  (d) 

Figure  6.  The  maximum  pressure  at  the  pinch-off  location  at  different  resolutions  using 
different  threshold  breakup  approaches  as  in  Fig.  4. 


The  threshold  breakup  tests  discussed  above  verified  that  the  reason  for  simulation 
divergence  is  the  singularity  present  at  the  pinch-off  instant.  Bypassing  the  singularity  can 
effectively  provide  a  grid-converged  simulation.  However,  the  remaining  question  is  how  to 
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develop  a  mathematically/physically  consistent  model  in  order  to  bypass  the  pinch-off 
singularity  and  still  accurately  capture  the  process  post  pinch-off.  This  motivates  the 
development  of  a  meso-scale  approach  to  achieve  a  detailed  understanding  of  the  pinch-off 
physics. 


1.2  Development  of  a  singularity-free  meso-scale  simulation 
approach 


As  shown  in  the  previous  section,  numerical  convergence  near  liquid  pinch-off  is  difficult  to 
achieve  due  to  the  singular  nature  of  the  continuum  formulation.  Physically,  the  effects  of 
molecular  scale  fluctuations,  which  are  averaged  out  in  the  continuum  formulation,  start  to  play 
a  critical  role  near  pinch-off.  To  correctly  account  for  all  the  physics,  a  new  model  needs  to  be 
developed  to  incorporate  molecular-scale  information.  Doing  so  in  the  continuum  framework  is 
challenging  because  the  framework  is  established  on  the  Navier-Stokes  equations  based 
continuum  assumption.  Particle  based  approach  such  as  molecular  dynamics  can  directly 
capture  the  molecular  events,  but  at  a  prohibitively  high  cost.  In  this  project,  the  mesoscopic 
Dissipative  Particle  Dynamics  (DPD)  approach  is  adopted  to  bridge  the  molecular  and 
continuum  scales.  The  interfacial  dynamics  is  accounted  for  using  a  variant  of  DPD  called 
Many-body  DPD  described  below. 


1.2.1  Background  of  MDPD 

The  advantage  of  dissipative  particle  dynamics  (DPD)  resides  in  the  simplicity  of  the  underlying 
algorithm  of  particle  interaction  under  a  soft  repulsive  potential.  It  has  been  shown  that  DPD  can 
be  constructed  as  a  mesoscopic  model  of  molecular  dynamics  (MD)  [6],  where  the  coarse- 
graining  of  Lennard- Jones  clusters  can  lead  to  a  suitable  DPD  force  field.  This  result  holds  at 
sufficient  low  densities  and  for  small  ranges  of  gyration,  when  many-body  effects  can  be 
neglected. 

DPD  has  been  used  to  investigate  phase  separation  in  immiscible  binary  liquid  mixtures  [7] -[9] 
droplet  deformation  and  rupture  in  shear  flow  [10],  and  droplets  on  surfaces  under  the  influence 
of  shear  flow  [11].  In  single-species  fluid  problems,  however,  the  standard  DPD  method 
presents  a  fundamental  limitation,  in  that  the  repulsive  soft  potential  alone  cannot  reproduce 
surface  tension.  This  potential  leads  to  a  predominantly  quadratic  pressure-density  equation  of 
state  (EOS)  [12],  while  a  higher-order  pressure-density  curve  is  necessary  for  the  coexistence 
of  the  liquid  and  vapor  phases. 

The  free  energy  of  the  DPD  system  was  modified  by  Tiwary  and  Abraham  [13]  to  include  a 
correction  directly  derived  from  the  van  der  Waals  EOS.  The  new  term  added  the  long-range 
attractive  force  that  is  responsible  for  surface  tension  between  the  liquid  and  the  vapor  phase; 
surface  tension  emerges  from  the  asymmetry  of  the  intermolecular  forces  acting  on  a  layer  of 
molecules  at  the  liquid-vapor  interface.  As  this  asymmetry  causes  larger  intermolecular 
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distances  in  the  outer  layer  than  in  the  liquid  bulk,  the  forces  in  the  layer  act  to  contract  the 
interface. 

The  many-body  DPD  (MDPD)  method  by  Pagonabar-raga  and  Frenkel  [14]  also  included  an 
attractive  force.  The  amplitude  of  the  soft  repulsion  was  made  proportional  to  the  local  density 
of  the  particles,  thus  achieving  a  cubic  pressure-density  relation.  A  similar  approach  was 
introduced  by  Nugent  and  Posch  in  the  context  of  smoothed  particle  hydrodynamics  (SPH)  [15]. 
The  connection  between  MDPD  and  SPH  was  then  clarified  by  Espanol  and  Revenga  [16],  who 
introduced  a  smoothed  DPD  (SDPD)  method  as  a  SPH  variant  based  on  a  new  formalism 
developed  for  discrete  hydrodynamics  [17]. 

MDPD  was  extensively  investigated  by  Warren  [18]  and  Trofimov  et  al.  [19].  The  method  can 
be  used  for  the  study  of  single  species  free-surface  flow,  for  instance,  in  the  case  of  pinch-off  of 
a  liquid  thread  during  the  formation  of  a  drop.  There,  the  vanishing  liquid  thread  diameter 
induces  a  singularity  in  the  hydrodynamic  description  of  a  two-phase  flow  interface  that  is  at 
odds  with  the  continuum  description  of  surface  tension  as  a  taut  membrane  of  zero  thickness. 
However,  the  focus  of  past  research  has  been  mostly  on  bulk  properties,  namely,  the  pressure- 
density  relation  [18],  This  project  focuses  instead  on  the  interfacial  behavior  of  the  MDPD  fluid, 
on  the  sensitivity  of  the  surface  tension  coefficient  a  to  the  particle  interaction  parameters,  and 
on  the  comparison  with  the  properties  of  real  liquids. 

In  the  following,  an  extensive  set  of  simulations  is  presented  for  a  liquid  with  number  density 
p  in  coexistence  with  a  very  dilute  vapor.  The  dependence  of  O  from  p  is  investigated  based  on 
a  table  of  representative  interaction  parameters.  An  approximate  mean-field  expression  is 
derived  to  directly  evaluate  the  surface  tension  coefficient  from  the  interaction  parameters. 
The  simulation  of  liquid  thread  pinch-off  is  then  verified  using  analysis  results  derived  in  the 
limit  of  a  slender  jet. 

All  the  simulations  are  enabled  by  the  particle  dynamics  software  code  LAMMPS  (Large- 
scale  Atomic/Molecular  Massively  Parallel  Simulator)  [20]  with  the  addition  of  a  new  MDPD 
class.  The  computationally  scalable  implementation  of  LAMMPS  guarantees  the 
optimization  of  the  particle  interaction  calculation  through  an  efficient  neighbor  list 
algorithm,  and  will  not  be  discussed  in  this  work. 


1.2.2  Numerical  details  of  MDPD 

MDPD  inherits  the  three  pairwise-additive  inter-particle  forces  formulation  of  the  standard 
DPD  scheme,  F.  =  ^  +  F(;  +  At~1/2F.f  ,  with  At  the  simulation  time  step.  The 

conservative,  dissipative,  and  random  forces  of  this  expression  are  defined,  respectively,  as 

K  =-r^(nMj-rij)ru, 
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where  r  =  r  lr.  and  v.  =  v,  -y , .  Warren’s  approach  [18]  is  pursued  for  F,f  ;  the  repulsive 

ij  ij  ij  lJ  1  J  x  x  A  y  A 

force  depends  on  a  weighted  average  of  the  local  density,  whereas  the  attractive  force  is 
density-independent, 


FS  =  A,wc  (r,j)  +  +  Pj  )wd  (r,) . 

The  weight  functions  wc(r)  =  (1  —  r/rc)  and  wd(r)  =  (1  —  r/rd)  vanish  for  r  >  rc  and  r  >  rd, 
respectively.  Since  a  DPD  method  with  a  single  range  may  not  have  a  stable  interface,  the 
repulsive  contribution  is  set  to  act  at  a  shorter  range  rd  <  rc  than  the  soft  pair  attractive 
potential. 

The  many-body  repulsion  is  chosen  in  the  form  of  a  self-energy  per  particle  which  is 
quadratic  in  the  local  density,  +  Pj)wd (rtj) ,  where  B  >  0.  The  density  for  each  particle  is 

defined  as 


Pi  =2>pA,)’ 

j*i 


and  its  weight  function  vvp  is  defined  as 

wp{r)=^j(l-rlrd)2 

Inr 

wp  vanishes  for  r>  rd  and  for  convenience  is  normalized  so  that  J d3rwp  (r)  =  1 . 

The  DPD  thermostat  consists  of  random  and  dissipative  forces.  The  6 '{j  coefficients  are 

independent  identically  distributed  Gaussian  random  numbers  with  zero  mean  and  unit 
variance.  The  equilibrium  temperature  T  is  maintained  through  the  condition  posed  by  the 
fluctuation-dissipation  theorem 


?=2)kBT, 


where  kg  is  the  Boltzmann  constant.  The  weight  function  for  the  dissipative  force  is 

wD(r)  =  (l-r/rc)2. 

All  the  simulations  presented  in  this  report  are  carried  out  with  the  velocity  Verlet  algorithm 
of  Groot  and  Warren  [21]  using  the  value  0.5  for  the  empirical  parameter. 

Dissipative  particle  dynamics  methods  operate  in  reduced  units,  such  that  energy  is  measured  in 
units  of  kgT\  length  in  units  of  rc,  and  mass  in  units  of  mass  of  a  single  particle,  m.  In  the 
following  discussion,  we  will  find  it  convenient  to  refer  to  dimensionless  quantities  such  as  the 
Ohnesorge  number  and  the  Bond  number. 
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1.2.3  Application  to  liquid-vapor  interface  dynamics 


The  operational  expression  to  calculate  surface  tension  at  liquid-vapor  equilibrium  uses  the 
macroscopic  normal  and  tangential  pressure  values,  or,  alternatively,  a  perturbation  formalism 
where  surface  tension  is  expressed  as  the  energy  difference  between  a  reference  state  and  a 
state  with  an  infinitesimal  area  variation.  Reference  [21]  compares  the  thermodynamic  and 
mechanical  routes  for  MDPD.  We  examine  here  the  simplified  relation  that  is  established 
between  the  pressure  difference  AP  across  the  liquid-vapor  interface  and  the  shape  of  a 
capillary  surface  according  to  the  Young-Laplace  equation  (Y-L), 


AP  =  cr(—  +  — ) 
R,  R2 


In  this  expression,  R/  and  R2  are  the  principal  radii  of  curvature  of  the  surface.  The  Y-L  relation 
is  deduced  from  mechanical  stability,  and  it  has  been  shown  to  hold  even  for  nanometer-size 
bubbles  [23].  The  internal  pressure  of  a  liquid  can  be  calculated  from 

where  (  )  denotes  the  sample  average  of  the  particles  contained  in  the  volume  V  and  v  is  the 

sample  velocity  after  being  spatially  averaged  on  V.  The  first  term  on  the  right-end  side 
represents  the  thermal  agitation  of  the  molecules  of  the  system.  The  second  term  is  due  to  the 
interaction  potential. 

By  way  of  illustration,  drops  with  different  diameters  are  simulated  from  an  initial  spherical 
lattice  of  particles  in  a  box  with  periodic  boundary  conditions.  The  size  of  the  box  ensures  that 
the  drop  is  isolated.  The  initial  random  velocity  distribution  is  at  kBT  =  1.  This  value  is 
maintained  in  the  simulation  by  fixing  the  parameters  £  =  1  and  y  =  0.5.  The  virial  pressure  at 
equilibrium  is  plotted  in  Fig.  7  as  a  function  of  radius  for  the  set  A  =  -40,  B  =  25,  rc  =  1,  rci  = 
0.80.  In  this  exercise  there  is  no  change  in  coarse-graining  when  the  droplet  radius  is  varied 
from  one  simulation  to  the  next. 

In  examining  the  diagram  in  Fig.  7(a),  it  is  apparent  that  the  thermal  contribution  of  the  virial 
pressure  is  constant  over  most  of  the  drop  and  that  the  plateau  value  is  pkBT  (the  number 
density  is  4.76).  At  the  liquid-vapor  surface,  the  thermal  contribution  decreases  to  zero  within 
approximately  one  unit  length.  It  can  be  shown  that  the  actual  extent  of  this  transition  region 
depends  on  the  temperature  of  the  drop  and  becomes  larger  at  higher  temperatures.  The 
conservative  term  in  Fig.  7(b)  has  a  more  complex  behavior.  Past  the  fluctuations  near  the 
droplet  center,  this  term  reaches  a  constant  value  when  the  number  of  sampled  particles 
becomes  sufficiently  large. 

The  negative  dip  at  the  drop  surface  can  be  attributed  to  the  strengthening  of  the  attractive  force 
as  outward  layer  particles  on  average  possess  a  larger  inter-particle  distance  compared  to  the 
bulk  particles.  Further  away,  the  longer-range  attractive  force  vanishes,  the  virial  pressure 
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becomes  less  negative,  and  it  eventually  goes  to  zero  when  the  average  inter-particle  distance  is 
larger  than  rc. 

The  sum  of  the  conservative  and  thermal  terms  forms  the  overall  pressure  field  and  is  displayed 
in  Fig.  7(c).  The  difference  between  the  normal  and  the  tangential  component  of  the  spherical 
interface  is  displayed  in  Fig.  7(d).  For  clarity,  the  diagrams  from  the  smallest  of  the  three 
droplets  are  omitted  in  this  figure.  The  positive  peak  at  the  drop  interface,  followed  by  noise  in 
the  bulk  (also  omitted  for  clarity),  illustrates  the  non-hydrostatic  nature  of  the  virial  stress 
tensor  that  is  responsible  for  surface  tension.  This  diagram  is  consistent  to  what  has  been 
observed  in  molecular  dynamics  simulations  for  a  Lennard- Jones  fluid  in  a  slab  geometry  [24], 

The  surface  tension  coefficient  can  be  calculated  directly  according  to  the  Y-L  relation  from  the 
slope  of  the  lines  fitting  the  values  of  AP  at  various  diameters;  see  Fig.  8.  It  is  noted  that  this  is 
not  the  only  methodology  to  evaluate  a  and  that  a  more  complete  discussion  can  be  found,  for 
instance,  in  Refs.  [13]  and  [21].  The  radius  of  curvature  of  the  drop  is  determined  from  the 
point  where  the  number  density  falls  below  half  of  the  bulk  density  value.  This  is  a  relatively 
precise  measurement  given  the  steepness  of  the  interface  at  this  temperature.  The  internal 
pressure  is  taken  from  the  virial  plateau  value,  with  an  uncertainty  of  +0.4  for  the  smallest 
droplets  and  of  less  than  +0.002  for  the  largest  drops  in  the  range  considered  here.  The  MDPD 
parameters  and  the  corresponding  values  of  a  are  listed  in  Table  I.  The  values  found  for  Sets  2 
and  7  are  in  close  agreement  with  the  values  reported  in  Ref.  [18]  (4.95  and  7.54)  for  the  same 
parameters. 


Set 

A 

B 

rd 

P 

0 

O  fit 

1 

-40 

40 

0.80 

3.94 

1.90 

3.33 

2 

-40 

40 

0.75 

5.12 

4.69 

5.63 

3 

-30 

25 

0.75 

5.27 

3.51 

4.48 

4 

-40 

80 

0.70 

5.47 

6.89 

6.36 

5 

-50 

40 

0.75 

5.60 

8.42 

8.43 

6 

-30 

20 

0.75 

5.80 

4.48 

5.44 

7 

-40 

25 

0.75 

6.10 

7.30 

8.02 

8 

-40 

40 

0.70 

6.51 

10.2 

9.11 

9 

-40 

20 

0.75 

6.77 

9.14 

9.90 

10 

-40 

25 

0.70 

7.82 

14.7 

13.2 

TABLE  I.  Parameter  sets  and  properties  of  MDPD  fluids.  The  bulk  density  p  and  the 
surface  tension  coefficient  <7  are  calculated  for  kuT  —  1. 
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(a) 


Radius 


(b) 


(c)  (d) 

Figure  7.  Virial  pressure  as  a  function  of  radius  for  three  drops  of  fluid  1  (the 
interaction  parameters  are  listed  in  Table  I):  (a)  contribution  of  thermal  motion;  (b) 
contribution  of  the  conservative  forces  (c)  resulting  virial  pressure;  (d)  radial  minus 
tangential  component  of  virial  pressure. 
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Figure  8.  Pressure  difference  across  the  liquid-vapor  interface  of  a  drop  in  equilibrium 
versus  its  curvature  for  the  10  MDPD  fluids  listed  in  Table  I. 


The  occurrence  of  capillary  pinch-off  has  already  been  demonstrated  by  a  modified  DPD 
method  [24]  as  well  as  by  pioneering  molecular  dynamics  simulations  [26],  Here,  we 
examine  whether  MDPD  captures  the  correct  fluid  dynamic  behavior  by  assessing  the 
response  of  a  quiescent  circular  jet  of  diameter  D  to  a  periodic  axial  perturbation.  It  is  well 
known  that  the  jet  becomes  unstable  if  A/D  >  n  [27].  Above  that  cutoff  wavelength,  the 
cylinder  surface  is  driven  toward  the  state  with  the  smallest  area,  and,  locally,  to  pinch-off.  In 
this  process,  the  greater  the  cohesion  between  the  particles  (and  thus  the  surface  tension),  the 
faster  the  instability  develops.  Below  the  cutoff,  surface  tension  has  the  opposite  role  of 
restoring  the  jet  to  its  unperturbed  state.  In  the  inviscid  description  of  the  Plateau-Rayleigh 
instability  [28],  linear  analysis  of  perturbation  leads  to  the  maximum  growth  of  a  disturbance 
for  A/D-4.5.  By  dropping  viscosity  from  dimensional  analysis,  the  characteristic  pinch-off 
time  is 

T  =  ^~^Ja 

As  a  demonstration,  simulations  are  carried  out  for  two  jet  diameters  in  a  103  box  with 
periodic  boundary  conditions  and  sinusoidal  perturbation  of  wavelength  A  =  10.  The  MDPD 
parameters  correspond  to  Set  2  in  Table  I.  The  random  and  dissipation  coefficients  are  t,  =  12 
and  y  =  72;  the  time  step  is  At  =  0.001. 

In  the  first  simulation,  we  set  L/D  =  4.5  with  251  particles.  Fig.  9(a)  displays  a  sequence  of 
snapshots  as  seen  through  a  slice  of  the  jet.  Initially,  the  amplitude  of  the  perturbation  is  only 
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one  layer  of  particles.  The  jet  begins  to  neck  at  around  t  =  2.7  and  pinches  off  soon  thereafter, 
so  that  by  t  =  13  the  two  half  droplets  have  fully  equilibrated.  The  characteristic  time  is  x  = 
3.46.  The  actual  time  to  pinch-off  depends  on  the  seed  values  that  are  used  to  generate 
random  numbers  for  the  stochastic  force  and  for  the  initial  velocity  distribution.  It  is, 
however,  well  defined  in  every  realization  of  this  case. 

In  the  second  simulation,  the  diameter  is  increased  to  D  =  5,  with  1116  particles.  No  capillary 
instability  develops,  even  for  large  amplitude  perturbations,  and  the  jet  returns  to  a  cylindrical 
shape  under  the  restoring  force  of  surface  tension.  The  simulation  can  continue  for  a  very 
long  time  (the  last  frame  of  Fig.  9(b)  shows  a  snapshot  at  t  =  200)  without  pinch-off  ever 
occurring.  This  and  the  previous  result  are  consistent  with  the  Plateau-Rayleigh  analysis. 


Figure  9.  Snapshots  of  capillary  dynamics  of  a  liquid  thread:  L/D  =  4.51  (top)  and  L/D  = 
2.5  (bottom). 


1.2.4  Non-periodic  boundary  conditions 

A  novel  element  in  this  study  is  the  implementation  of  a  Non-Periodic  Boundary  Condition 
(NPBC)  to  broaden  the  application  of  particle  methods  and  study  straining  gas-liquid  fields. 
A  first  challenge  is  that  the  boundary  conditions  need  to  include  the  right  amount  of  attractive 
force,  responsible  for  the  continuity  of  the  liquid.  Implementation  of  a  simple  reflecting  wall 
alone  is  not  sufficient  because  there  would  not  be  any  force  anchoring  the  MDPD  boundary 
particles.  Also,  an  effective-force  approach  where  potential  is  applied  to  impose  the  desired 
kinematics  is  insufficient  to  keep  the  liquid  and  gas  phase  separated. 

As  shown  in  the  schematics  of  Fig.  10,  two  layers  of  particles  are  built  into  the  domain  on 
each  side  of  the  computational  box  where  the  NPBC  is  assigned.  The  outermost  layer  (O)  is 
modified  at  every  iteration  by  placing  particles  of  the  prescribed  type.  This  layer  is  a  fixed 
(no  time  integration)  barrier  whose  composition  depends  on  the  instantaneous  location  of  the 
boundary  interface.  Layer  (O)  provides  the  necessary  attractive  force  for  the  continuity  of  the 
MDPD  liquid. 
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Particles  from  the  innermost  (I)  layer  of  the  buffer  are  free  to  move  according  to  the 
distribution  of  the  surrounding  particles.  At  the  end  of  every  iteration  their  velocity  is 
modified  so  that  1)  the  average  particle  velocity  of  a  cell  of  the  buffer  matches  the  prescribed 
velocity  of  that  buffer;  2)  the  squared  deviation  from  that  average  value  matches  the 
prescribed  temperature.  Here,  “cell”  refers  to  an  element  of  volume  of  side  U=  nf  ,  where  n 
is  the  number  density  of  component  i;  see  Fig.  10.  Particles  are  free  to  cross  in  both  directions 
the  boundary  between  the  (I)  layer  and  the  domain  proper,  but  particles  crossing  the 
separation  between  layers  (I)  and  (O)  are  reflected  back.  This  wall  moves  at  the  velocity  of 
the  normal  component  of  v*. 


side  normal 


inner  layer 

%  $--* 

*o  Id1- 

outer  layer 


P 


Figure  10.  Schematics  showing  the  prescription  of  composition  and  velocity  at  a  face 
boundary. 

A  numerical  test  discussed  here  is  the  enforcement  of  the  extensional  configuration  displayed 
in  Figure  2,  with  minimum  radius  of  the  dumbbell  R  =  6.5.  The  flow  field  is  achieved  by 
imposing  NPBCs  on  the  sides  of  a  503  box.  The  domain  initially  contains  14,512  particles  of 
liquid  and  15,924  particles  of  gas.  The  liquid  particles  have  the  properties  listed  above;  the 
gas  particles  are  treated  as  DPD  particles  (no  attractive  force):  A  =  10,  rc  =  1,  B  =  0,  rd  = 
0.75;  y=  288,  and  £=  24.  These  parameters  correspond  to  v=  0.84  and  n  =  1.2.  To  obtain  a 
liquid/gas  density  ratio  of  780,  the  mass  is  set  to  m  =  0.0064. 

The  time-averaged  velocity  of  the  particles  arranges  itself  quickly  to  a  fully  developed 
extensional  field.  As  time  progresses,  liquid  particles  exit  the  domain  while  more  gas 
particles  enter  it.  To  account  for  the  advection  component  involved  in  imposing  a  prescribed 
strain  rate,  the  time  step  is  set  to  0.00004  in  DPD  units.  Thus,  a  calculation  including  the 
interaction  with  the  gas  phase  is  substantially  more  expensive  than  the  single-component 
simulations  of  capillary  pinch-off  discussed  earlier.  Also,  NPBCs  require  two  additional 
layers  of  particles  for  each  boundary  pair  -  a  12%  increase  in  the  overall  system  size.  Fig.  11 
shows  an  excellent  agreement  with  the  corresponding  CLSVOF  simulation  at  three  instants 
preceding  pinch-off. 
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Figure  11.  Axis-symmetric  dumbbell  subject  to  extensional  flow:  CLSVOF  (top)  and 
MDPD  (bottom)  snapshots  before  pinch-off. 
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1.3  Scaling  analysis  of  pinch-off  using  meso-scale  simulations 


1.3.1  Theoretical  background  of  pinch-off  scaling 

Theoretical  description  of  interface  dynamics  when  a  liquid  thread  shrinks  to  zero  near  pinch- 
off  -  and  interface  curvature  goes  to  infinity  -  has  been  rare  in  the  literature.  A  notable 
exception  is  Eggers’  similarity  solution  derived  for  drop  pinch-off  in  free-surface  flow  [29].  It 
is  shown  that  pinch-off  at  (xo,  to)  is  a  singularity  of  the  Navier-Stokes  equations  in  a  critical 
region  Ix'l  «1  and  t'«  1 ,  where 


x  = 


and 


t  = 


t-tn 


The  parameters  £  v  and  tv  depend  only  on  the  property  of  the  liquid, 

£v=(p\?)/cr  and  tv  =(p2v3)/a  , 

with  v,  p,  and  a  are  the  liquid  viscosity,  density  and  surface  tension,  respectively.  For  water 
(p  =  1000  kg/m3;  v=  1.138  10-6  m2/s;  a  =  0.0728  N/m),  these  viscous  length  (£v  =  1.77  10'8 

m)  and  time  (rr  =  2.78  10~10  s)  scales  are  very  small  compared  to  scales  in  most  continuum 

simulations.  An  important  consequence  of  Eggers’  work  is  that  the  minimum  liquid  thickness 
along  the  ligament  scales  linearly  with  the  time  left  to  the  singularity  time  to, 

hm,=  0.03—  (»-»„), 


Figure  12.  Sketch  of  the  flow  geometry  investigated  by  Eggers  [29]. 


Another  analytical  work  by  Lister  and  Stone  [3]  applied  the  asymptotic  balance  between 
surface  tension  and  viscous  stress. 


V/Kn  ~  PUV/hnin  , 

which  also  suggest  the  thinning  rate  of  the  thread  u  becomes  asymptotically  constant  in  time, 

u  -  (7/{pv) 
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and  pinching-off  occurs  with  a  finite  period.  All  the  work  suggests  that  Navier-Stokes 
equation  is  singular  at  the  point  of  pinch-off  and  it  has  a  self-similar  solution  in  the  limit  of 
the  singularity.  The  solution  is  universal,  in  the  sense  that  it  is  independent  from  the 
particular  initial  or  boundary  conditions  of  the  experiment. 

The  physics  of  liquid  pinch-off  is  in  fact  more  complex  than  that  described  by  Navier-Stokes 
equations.  Intuitively,  it  can  be  expected  that  when  the  liquid  thread  thickness  is  comparable 
to  the  scale  of  thermal  fluctuations  that  are  averaged  out  in  the  continuum  description  -  or  to 
the  actual  interface  width  (~100nm)  -  the  Navier-Stokes  description  of  two-phase  flow 
breaks  down.  For  liquid,  the  thermal  length  scale, 

£T=(kBT/c x)1/2, 

accounts  for  the  fact  that  molecules  at  the  interface  have  energies  due  to  both  thermal  motion 
and  surface  tension.  iT  is  of  the  order  of  a  nanometer  in  most  fluids  at  room  conditions.  In 
free  surface  flow,  a  threshold  length  exist, 


^  thres  ^  v 


UTV'5 

\^v  J 


below  which  thermal  capillary  instabilities  can  cause  the  interface  to  behave  stochastically. 
For  fuel  at  ambient  temperature,  £thres  =  1  pm,  which  is  comparable  to  the  thickness  of  the 
liquid  thread  that  forms  before  pinch-off. 

The  domain  concerned  by  thermal  fluctuations  (up  to  a  few  hundreds  of  nanometers)  has 
been  recently  probed  by  Molecular  Dynamics  (MD)  simulations.  MD  simulations  of  nanojets 
show  a  substantial  change  in  behavior  near  pinch-off  [30],  This  behavior  is  predicted  by  the 
stochastic  lubrication  equation  (SLE)  [26],  derived  from  the  Navier-Stokes  equations  in  the 
limit  of  a  slender  jet  and  with  the  addition  of  Gaussian  noise  to  the  deterministic  stress  tensor. 
The  relation  between  hmi„  and  to  -  t  is  described  by  a  power  law  with  exponent  0.418,  found 
by  numerical  integration  [31]. 


1.3.2  MDPD  scaling  analysis  for  capillary  pinch-off 

To  establish  a  connection  with  the  theory,  the  results  from  previous  MDPD  simulation  of 
dumbbell  drop  breakup  by  capillary  effects  (subsection  1.2.3)  are  post-processed  and  the 
minimum  thread  radius  hmin  are  plotted  as  a  function  of  the  time  to  breakup,  to  - t  in  Fig.  13. 
The  time  to  pinch-off  in  Fig.  13  is  normalized  by  x,  whereas  is  normalized  by  the 
unperturbed  jet  radius  R.  To  track  the  minimum  jet  radius,  it  is  necessary  to  post-process 
several  snapshots  of  particle  positions,  each  axially  sliced  into  50  bins.  It  is  assumed  that  the 
only  1%  or  less  of  the  particles  lie  outside  the  surface  of  the  jet,  a  threshold  consistent  with 
almost  no  vapor  phase.  The  jet  radius,  calculated  with  respect  to  the  instantaneous  position  of 
the  center  of  mass  in  each  slice,  is  defined  based  on  the  radial  histogram  of  number  density. 
The  pinch-off  time  is  established  as  the  instant  when  one  of  the  bins  becomes  empty.  The 
simulation  is  carried  out  in  a  803  periodic  cube  using  the  parameter  set  7  from  Table  I  for  a 
system  of  118,657  particles;  the  time  step  is  At  =  0.001.  The  dissipation  and  random 
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coefficients  are  a  =  0.5  and  t,  =  1,  giving  v  =  0.072;  the  Ohnesorge  number  is  therefore  Oh  = 

0.022. 


The  logarithmic  plot  of  hmin(to  ~  t)  in  Fig.  13  demonstrates  that  a  single  MDPD  simulation 
can  span  the  three  scaling  behaviors  listed  above.  Initially,  and  for  almost  a  time  decade,  the 
tracked  points  are  aligned  along  the  2/3  slope.  The  inertial-viscous  behavior  appears  at 
approximately  hmi„=0.3R  and  is  quickly  overtaken  by  a  trend  where  most  of  the  points  align 
along  the  slope  0.418.  This  occurs  approximately  at  hmin=0.15R,  above  the  thermal  capillary 
length  lT  =  0.042  R.  At  that  point,  as  noted  by  Eggers  [31],  any  random  fluctuation  which 
increases  the  thread  radius  will  also  increase  its  effective  mass,  slowing  down  the  motion  of 
the  fluid;  conversely,  any  fluctuation  toward  a  smaller  neck  radius  will  accelerate  pinch-off. 
The  gap  between  deterministic  hydrodynamics  and  molecular  dynamics  is  thus  bridged  in 
this  simple  example.  Other  MDPD  simulations  (not  shown  here)  confirm  that  a  larger 
viscosity  or  a  smaller  domain  causes  the  potential-flow  scaling  to  disappear  from  the 
diagram. 


Figure  13.  Variation  of  minimum  thread  radius  h„,in  versus  time  to  breakup  to  -  t.  The 
three  lines  correspond  to  potential  flow  (slope  2/3),  inertial -viscous  (slope  1),  and 
stochastic  (slope  0.418)  power  laws.  The  open  symbols  are  measured  values  from  the 
MDPD  simulation. 


1.3.3  MDPD  scaling  analysis  for  pinch-off  in  straining  flow 

To  simulate  pinch-off  of  drop  in  a  straining  flow  using  MDPD,  Non-Periodic  Boundary 
Conditions  (NPBCs)  is  implemented  at  the  boundaries  of  the  simulation  domain.  As  time 
progresses,  liquid  particles  exit  the  domain  while  more  gas  particles  enter  it.  The  number 
density,  mass  and  viscosity  of  these  two  fluids  are  arranged  so  that  the  gas-liquid  density 
ratio  is  0.0013,  and  that  the  kinematic  viscosity  ratio  is  15.  These  two  values  are  in  the  range 
of  typical  air/liquid  ratios  for  air  at  ambient  conditions. 

The  straining  case  presented  in  subsection  1.1.2  was  simulated  using  MDPD.  Several 
snapshots  for  the  simulations  at  two  different  strain  rates  are  shown  in  Fig  14.  The  time  is 
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scaled  by  DPD  unit  time.  As  the  strain  rate  increases,  the  global-scale  elongation  of  drop 
becomes  faster.  However,  closer  look  at  the  local-scale  neck  region  (rectangular  box)  reveals 
that  the  local  pinch-off  process  does  not  depend  much  on  the  strain  rate.  At  the  same  time  t-to 
relative  to  pinch-off,  the  thread  shapes  are  quite  similar  for  the  two  strain  rate  cases  even 
though  the  larger-scale  drop  shapes  are  different. 

The  minimum  thread  radii  hmin  under  the  two  strain  rates  are  extracted  from  the  simulation 
and  plotted  as  a  function  of  time  in  Fig.  15.  The  two  cases  show  certain  degree  of  discrepancy 
due  to  the  molecular-scale  fluctuations  when  pinch-off  is  approached.  However,  hm(n  for  both 
cases  shows  the  same  transition  from  inertial  scaling  to  inertial  viscous  scaling  and  to 
stochastic  scaling.  The  point  of  transition  seems  agree  as  well  for  the  two  cases.  The  results 
suggest  that  strain  rate  does  not  affect  the  evolution  of  thread  neck  near  pinch-off.  As  pointed 
out  by  Eggers  [29],  the  thread  behavior  approaching  pinch-off  becomes  self-similar  and 
depends  only  on  intrinsic  length  and  time  scales,  independent  of  external  conditions  such  as 
strain  rate.  The  results  in  Fig.  15  agree  with  such  theoretical  arguments.  The  theoretical 
argument  also  suggests  that  the  transition  point  between  inertial  and  viscous  scaling  is  only 
dependent  on  the  fluid  properties  (such  as  density,  viscosity  and  surface  tension)  and 
independent  of  strain  rate.  Such  condition-independent  transition  point  separates  the  physics 
occurring  in  different  scales.  At  a  scale  above  the  transition  point,  strain  rate  affect  larger 
scale  two-phase  flow  (see  Fig.  14).  At  a  scale  below  the  transition  point,  the  physics  is 
independent  of  strain  rate  and  only  affected  by  molecular-level  fluctuations.  It  can  also  be 
inferred  that  because  of  the  scale  separation,  the  molecular  fluctuation  does  not  impact  larger 
scale  dynamics. 


t  25.2  30.2  43.2  45.2 

t-to  -18  -13  0  2 


t  13.6  18.6  31.6  33.6 

t-t0  -18  -13  0  2 


Figure  14.  Snapshots  of  dumbbell  drop  breakup  in  a  straining  flow  at  different  strain 
rate:  ^  =  0.4xlCfs'1  (top)  and  y  =  1 06  s"1  (bottom). 
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Figure  15.  Variation  of  minimum  thread  radius  hmi„  versus  time  to  breakup  to  ~  t.  for 
different  strain  rate:  y  =  0.4x1  Cf  s'1  (circle)  and  y  =  1 06  s'1  (dot). 
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1.4  Development  and  demonstration  of  continuum-scale  closure 
model 


1.4.1  Closure  model  development  from  scaling  analysis 

The  analysis  in  subsection  1.3.3  gives  rise  to  a  key  finding  that  allows  developing  a  physics- 
based  closure  model  for  continuum  simulations. 

The  simulation  based  on  the  continuum  formulation  can  resolve  the  thread  neck  down  to  the 
inertial  to  inertial  viscous  scaling  transition  point.  Because  the  process  below  the  transition 
point  is  independent  of  external  flow  conditions,  the  transition  point  can  be  used  as  a 
universal  threshold  for  certain  liquid-gas  pair  with  fixed  fluid  properties.  The  threshold  is 
physics-based  and  independent  of  grid-resolutions.  The  value  of  the  threshold  can  be 
obtained  using  meso-scale  analysis  such  as  MDPD  simulations  and  applied  in  continuum 
simulations  as  the  criterion  to  stop  simulating  the  process  below  the  transition  point  and 
introduce  a  numerically  forced  breakup.  Because  of  the  condition-independent  nature  of  the 
threshold,  same  threshold  for  numerical  breakup  can  be  applied  for  any  pinch-off  processes 
occurring  in  the  simulation  domain  irrespective  of  their  local  strain  rate.  The  scaling 
transition  point  thus  becomes  the  only  parameter  of  the  closure  model. 

In  next  subsection,  a  generic  implementation  of  such  closure  model  in  CLSVOF  framework 
will  be  described.  A  converged  pinch-off  simulation  when  the  closure  model  is  applied  will 
be  demonstrated. 


1.4.2  Generic  implementation  of  the  closure  model 

Since  a  physics-based  threshold  has  been  determined  as  the  scaling  transition  point,  the 
implementation  of  the  closure  model  involves  developing  a  generic  procedure  to  identify  and 
cut  the  liquid  in  thread  neck  region  that  falls  below  the  threshold. 

The  identification  procedure  is  motivated  by  the  algorithm  for  liquid  filament  identification 
proposed  by  Kim  and  Moin  [4],  but  with  substantial  modifications  to  adapt  to  the  current 
CLSVOF  and  AMR  framework.  The  algorithm  is  illustrated  in  Fig.  16.  First,  all  cells  are 
grouped  in  gas  phase  (17  <  0),  liquid  phase  (77>  0),  and  liquid-gas  interface  (77>  0  and 
directly  adjacent  to  77  <  0  ),  respectively.  Then  for  each  interface  cell  Xi,  its  unit  normal  iij  is 
calculated  as  n(  =  V ^/|V <f>\ .  With  Xi  and  iii,  a  traverse  line  can  be  constituted  that  serves  as 

the  center-line  of  a  three-dimensional  cylinder  traversing  through  the  liquid  phase  in  direction 
ni  with  a  radius  r.  Starting  from  the  interface  cell  Xj,  all  its  immediate  neighbors  in  the  liquid 
phase  are  searched.  If  the  location  of  a  neighboring  cell  falls  within  the  range  of  the  cylinder, 
it  is  added  into  a  list  of  cells.  The  next  round  search  initiates  from  the  immediate  neighbors  of 
the  newly  identified  cells  from  the  previous  round,  and  is  carried  out  in  the  way  of  a  multi¬ 
branch  tree  search  [5].  Note  that  the  cells  visited  in  one  round  of  search  are  tagged  such  that  a 
redundant  search  during  the  next  round  for  these  cells  can  be  avoided.  Eventually  the  search 
is  stopped  when  none  of  the  neighbors  of  the  cells  within  the  list  can  be  found  in  liquid  phase 
(including  liquid-gas  interface)  and  enveloped  by  the  cylinder.  Once  the  search  is  finished, 
the  maximum  distance  from  the  initial  interface  cell  x,  to  another  interface  cell  xe  in  the  list  is 
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calculated.  This  distance  represents  the  local  liquid  thread  thickness  in  the  direction  iij.  If  it 
does  not  exceed  the  threshold  pinch-off  size,  this  list  of  cells  is  marked  as  valid  candidate 
cells,  which  are  labeled  in  green  color  in  Fig.  16. 

The  above  cylinder-envelope  search  procedure  is  carried  out  for  each  interface  cell.  As  a 
result,  all  the  cells  meeting  the  criteria  of  liquid  thread  thickness  are  marked.  The  next  step  is 
to  collect  contiguous  cells  from  the  tagged  ones  into  separate  cell  lists  to  represent  separate 
liquid  structures.  This  is  achieved  again  by  performing  the  multi-branch  tree  neighbor  cell 
search  described  above.  Note  that  the  above  algorithm  is  generic  in  the  sense  that  any  region 
with  thickness  below  the  threshold  will  be  identified,  irrespective  the  location  and  orientation 
of  the  region.  The  algorithm  does  not  involve  numerical  computation  of  gradients,  making  it 
less  sensitive  to  numerical  perturbations.  Also  the  algorithm  can  be  applied  in  any 
dimensions. 

After  the  potential  liquid  pinch-off  structures  are  identified,  they  are  constrained  by  an 
additional  criterion:  A  liquid  filament  is  disqualified  if  it  contains  any  liquid  cells  that  touch 
the  boundaries  of  an  AMR  box.  In  other  words,  the  ligament  is  not  allowed  to  cross  an  AMR 
box  and  operations  on  the  ligament  are  constrained  within  one  local  box.  As  a  result,  the 
breakup  operation  may  be  delayed  in  case  the  ligament  does  cross  multiple  boxes.  This  might 
seem  arbitrary,  but  due  to  the  dynamic  features  of  AMR,  the  current  box-crossing  liquid 
structure  has  a  good  chance  to  be  completely  contained  by  a  new  box  within  a  few  time  steps 
as  the  mesh  refinement  needs  to  be  performed  adaptively  at  every  time  step.  In  addition,  it  is 
advantageous  to  consider  only  the  local  box  operation,  such  that  the  complexity  of 
exchanging  liquid  structure  information  between  boxes  computed  on  different  processors  is 
avoided,  thus  the  communication  overhead  is  saved. 

In  the  current  model,  the  forced  breakup  occurs  at  a  very  small  scale.  It  is  thus  reasonable  to 
assume  that  the  continuum  simulation  (CLSVOF  algorithm)  directly  captures  the  small  scale 
pinch-off  between  the  smallest  satellite  drops.  The  resultant  small  amount  of  residual  mass 
due  to  the  forced  breakup  is  considered  negligible;  albeit  the  mass  should  be  redistributed 
back  to  the  parent  drops  in  the  strict  sense  of  mass  conservation. 


<f>  <0 


Figure  16.  A  schematic  of  the  liquid  thread  identification  algorithm. 
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1.4.3  Demonstration  of  simulation  convergence  using  the  closure  model 

With  the  closure  model  embedded  in  the  continuum  solver,  a  three-dimensional  drop  breakup 
in  a  straining  flow  is  simulated  at  different  resolutions  (controlled  via  AMR).  Although  the 
flow  is  largely  two-dimensional  axisymmetric,  we  perform  the  simulation  in  the  3D  version 
of  the  solver  to  test  the  generality  of  the  closure  model  implementation.  The  cost  of  this 
demonstration  is  much  higher  than  the  test  cases  presented  in  subsection  1.1.2. 

The  physical  length  scale  of  the  scaling  transition  point  is  on  the  order  of  micro  meters.  To 
have  the  continuum  solver  capture  such  small  scales  within  affordable  simulation  time,  we 
hereby  perform  the  demonstration  for  the  pinch-off  of  a  small  dumbbell  drop  of  initial 
unperturbed  cylindrical  diameter  D  =8  pm,  wavelength  X  =  20  pm,  and  wave  magnitude  s  = 
2pm.  For  water-air  system,  the  physical  scaling  transition  point  is  found  to  be  1.72  pm  based 
on  MDPD  simulation  results.  The  continuum  simulations  are  performed  in  a  48  pmx48 
pmx48  pm  domain  with  coarsest  resolution  being  used  to  be  96x96x96.  The  coarsest  grid 
size  is  0.5  pm,  far  below  the  threshold.  That  means  the  drop  dynamics  before  the  threshold  is 
reached  can  be  resolved  by  the  simulation.  The  grid  setup  for  refined  cases  is  summarized 
below. 


Case  No. 

Base  grid 

Refinement  level 

Refined  grid  size 

1 

96x96x96 

0 

0.5  pm 

2 

64x64x64 

1 

0.375  pm 

3 

96x96x96 

1 

0.25  pm 

4 

64x64x64 

2 

0.188  pm 

The  snapshots  of  simulations  prior  to,  near  and  post  pinch-off  are  shown  in  Fig.  17  with  the 
results  from  different  resolutions  compared.  By  enforcing  breakup  at  the  same  physics-based 
threshold  independently  from  grid  resolution,  the  simulation  converges  with  decreasing  Ax. 
The  drop  shapes  match  quite  closely  in  the  three  refined  cases  with  AMR  (Ax  =  0.375  pm, 
0.25  pm,  0.188  pm).  The  post  pinch-off  shapes  match  each  other  at  the  same  time  instant.  To 
quantitatively  verify  the  convergence,  the  breakup  time  to  and  the  post  pinch-off  separation 
distance  dsp  are  extracted  from  the  simulations  and  plotted  as  a  function  of  grid  size  in  Fig. 
18.  It  is  observed  that  both  to  and  dsp  reaches  a  constant  value  as  the  grid  size  decreases.  The 
results  clearly  demonstrate  the  convergence  of  the  calculation,  which  can  finally  be 
considered  as  a  Direct  Numerical  Simulation  of  the  flow. 
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t  =  3.6 


t  =  4.0 


Figure  17.  Snapshots  of  drop  shape  predicted  by  simulations  at  different  resolutions 
using  the  physics-based  closure  model. 


Figure  18.  The  grid  convergence  of  pinch-off  time  (a)  and  drop  separation  distance  post 
pinch-off. 
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